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R-modes of a rotating neutron star are unstable because of the emission of gravitational radiation. 
We explore the saturation amplitudes of these modes determined by nonlinear mode-mode coupling. 
Modelling the star as incompressible allows the analytic computation of the coupling coefficients. 
All couplings up to n = 30 are obtained, and analytic values for the shear damping and mode 
normalization are presented. In a subsequent paper we perform numerical simulations of a large set 
of coupled modes. 
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I. INTRODUCTION 

In 1998, Andersson 0, and Friedman and Morsink 
|3| showed that r-modes of a rotating neutron star 
can be unstable via the gravitational-radiation driven 
Chandrasekhar-Fricdman-Schutz (CFS) mechanism Q 
In inviscid stars this mechanism causes any mode 
that is retrograde in the co-rotating frame but prograde 
in the inertial frame to grow as it emits gravitational 
radiation. 

For neutron star rotational frequencies ^ 100 — lOOOHz 
gravitational radiation from an unstable r-mode falls 
within the range in which ground-based gravitational 
wave detectors such as LIGO are sensitive. If detectable 
this signal would probe the internal structure of neutron 
stars. Moreover, it has been suggested that the CFS 
instability of the r-mode could limit neutron star rota- 
tional frequencies Q. The key question, which currently 
remains unanswered, is to what amplitude these modes 
grow and whether they can be detected at all. 

The slow growth rate of these secularly unstable 
modes, as well as their three-dimensional nature, makes 
the problem ill-suited to direct hydrodynamical simula- 
tion. Gravitational radiation slowly drains away the ro- 
tational energy of the star over a time scale well in ex- 
cess of 10^ rotation periods, too long a time period for 
numerically stable hydrodynamical simulations of a ro- 
tating neutron star. 

There have been three attempts at direct hydrodynam- 
ical calculations 0| These simulations (1) assume 
that the background star has uniform initial rotation, 
(2) use approximate linear eigenfunctions as initial per- 
turbations, and (3) assume either a large initial ampli- 
tude, or else artificially increase the radiation reaction 
force enough to ensure that a large amplitude is reached 
relatively quickly. 

In one direct hydrodynamical calculation, Stergioulas 
and Font used 116^^ grid points, and perturbed an 
equilibrium model with an excitation corresponding to an 
approximate linear I = m = 2 r-mode eigenfunction with 
a substantial amplitude (a — 1). They evolved the sys- 
tem in a fixed, static spacetime geometry for 26 r-mode 
periods, after which the amplitude decreased because of 



numerical viscosity. No energy leakage out of the r-mode 
was observed, so they concluded that nonlinear couplings 
would not inhibit the instability. 

Lindblom et al J\ did a Newtonian hydrodynamical 
simulation on a 64 x 128 x 128 grid, with the gravita- 
tional radiation added as a radiation reaction force. How- 
ever, they boosted the radiation reaction force by a fac- 
tor of 4500 to force the instability to grow substantially 
in the ~ 30 rotation periods spanned by their calcula- 
tion. The initial condition was a uniformly rotating star 
with a small amplitude perturbation proportional to the 
vector spherical harmonic Re{Y^). In their simulation, 
the mode amplitude grew exponentially at first, reach- 
ing an amplitude of order unity. After that, the growth 
was halted by dissipation resulting from breaking waves 
(shocks). 

Gressman et al Q performed a Newtonian hydrody- 
namical simulation for non- viscous fluid flow very similar 
to 01 using a high resolution shock capturing scheme to 
accurately integrate their equations. They began by am- 
plifying an initially weak {I — m — 2) r-mode at 9 x 10^ 
times the radiation reaction rate, to amplitudes of or- 
der unity, at which point they ceased driving the mode 
and observed its decay. Subsequently the r-mode de- 
creased in amplitude slowly for a short time, but then 
decayed catastrophically. Gressman et al interpreted the 
dropoff as a slow leakage of energy to other fluid modes, 
which then act nonlinearly with the r-mode, resulting 
in the rapid decay. The time to decay was found to 
be amplitude-dependent: the larger the initial amplitude 
the faster the eventual decay. Previously, we have argued 
that this behavior can be attributed to the non-linear pe- 
riod of a particular three- mode coupling 0. 

None of these simulations excludes the possibility that 
the r-mode saturates at low amplitudes. If the timescale 
for the growth of the (I = m — 2) r-mode to saturate 
by nonlinear hydrodynamical effects is longer than the 
unphysically shortened growth rates used in the simula- 
tions, saturation at small amplitudes could be prevented 
artificially. Moreover, the spatial resolution of the hydro- 
dynamic codes may not be sufficient to resolve the small 
scale modes to which the r-mode couples effectively, ar- 
tificially increasing the possible r-mode amplitude. 
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To gain a deeper understanding of the r-mode instabil- 
ity and its saturation amplitude, we investigate the non- 
linear coupling of a large network of the inertial modes 
of a star. As long as the amplitudes of the modes remain 
small, it is a good approximation to treat the oscillations 
with (weakly nonlinear) perturbation theory. By treat- 
ing the neutron star as a uniform density incompressible 
object, the modes and their couplings can be calculated 
analytically. Our dynamical model then reduces to the 
simultaneous solution of a network of equations of the 
form 10] 

CA{t) - iWACA + lACA = -i— ^ K-^g^^CBCc (1) 

BC 



Arras et al based their WKB analysis on the Bryan 
modes to determine the scalings of the damping, normal- 
izations and couplings in their weak turbulence picture of 
r-mode saturation. Whereas they employed the Cowling 
approximation in computating the coupling coefficients, 
we do not. In fact, the perturbations of the gravitational 
potential make a nonnegligible contribution to the cou- 
pling coefhcients in this problem. 

Arras et al also presented a picture where a cascade 
of energy from the r-mode to smaller length scales is set 
up, and an inertial range is established. Our papers are 
designed to test this picture "experimentally" , and to see 
whether an alternative scenario is valid where the r-mode 
is effectively damped by coupling to just a few modes. 



Here the amplitudes \ca\ are assumed to be small, so 
we retain only the nonlinear couplings among triplets of 
modes, which, in the absence of dissipation and driving, 
would form an infinite dimensional Hamiltonian system 
with three- mode couplings k-jbc- Gravitational radi- 
ation and viscosity are added as driving and damping 
terms 7^. The quantities wa and sa are the mode fre- 
quency and mode energy of mode A respectively and a 
bar denotes complex conjugation. 

This approach has the great advantage that we can 
track the energy transfer and identify the couplings that 
become important. However it also poses many analyt- 
ical and computational challenges. We present the ana- 
lytic foundation to this investigation in this paper. The 
numerical investigation of the evolution of the oscillator 
network and the saturation amplitudes achieved will be 
presented in a later paper pjl |. 

First the linear eigenfunctions have to be obtained and 
appropriately normalized. In the case of a uniformly ro- 
tating uniform-density star an elegant coordinate trans- 
formation introduced by Bryan yields an exact analytic 
solution 10 A few properties of these generalized 

r-modes, their frequency structure and analytic normal- 
ization are reviewed in Section Hll 

Section IIIII details damping and driving terms of the 
linear modes. The analysis is a brief summary of that 
presented by Lockitch and Friedman 14] and Lindblom, 
Owen and Morsink and serves as a partial test of our 
approach to computing integrals over the star. Analytic 
scalings for the damping are obtained in the zero rotation 
limit. 

The main body of this paper explores a method of com- 
puting the nonlinear mode- mode couplings. In Section 
l iyi we specialize the formalism presented by Schenk et al 
|lQ,] to the case of the inertial modes of an incompress- 
ible star. Computing the coupling coefficients for modes 
with large mode number n = I + 1 is tricky because of 
precision problems. A technique for computing these co- 
efficients to arbitrary precision is presented in Section 
Ivl Computational time however limits us to computing 
only coefficients up to n = 30. Section IvTl explores some 
of the properties of the coupling coefficients, including a 
complete list of selection rules they obey. 



II. EQUILIBRIUM SOLUTION AND LINEAR 
MODES 



Consider a rotating ball of liquid with constant density 
that obeys the Euler equation 



dtv + V -Vv = -V I - - $ 



the continuity equation 

dtp + V-ipv) =0 
and Newtonian gravity 

= -iirGp 



(2) 



(3) 



(4) 



Here, v is the velocity field, p the density, $ the gravi- 
tational potential of the star and p the pressure. 

Maclaurin spheriods are uniform density, uniformly 
rotating solutions with eccentricity e and volume V = 
^ttR^. The average radius of the spheroid R is related to 
the equatorial (Re) and polar (Rp) radii of the spheroid 
by R^ = R^Rp. The pressure and angular velocity of the 
spheroid can be expressed in terms of the eccentricity as 
follows: 



= 2ttGp I -^^ ^ arcsme ^{1 - e'') I 



(6) 

where G is Newton's Constant and — x"^ + y'^ is the 
cylindrical radial coordinate. For a given eccentricity, Pe, 
Re and Rp are constant: 



Pe = 2ttGp' 



Re — R 



1 



1/6 



1 -e2 



Rp = i?(l-e2) 



■ arcsinej (7) 

(8) 
(9) 
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These solutions are dynamically and secularly sta- 
ble below the bifurcation point aX e — 0.81267, Q. — 
0.61174v^i^ XI}. 



Here a is an arbitary constant and Co = (1 — e^)/e^. The 
subscripts s^o and sfo indicate the exterior and interior 
regions of the star respectively. The coordinate transfor- 
mation to spheroidal coordinates is 



A. Linear Eigenmodes 

The linear stability analysis was carried out by Bryan 
in 1889 |l2|. Ipser and Lindblom |0 have given a mod- 
ern and improved rendering of this analysis and it is 
mainly their approach and notation on which we rely. 

The linearized perturbation equations for the La- 
grangian displacement ^ are 



UJ 



i + 2nxi^ -\/su 

V2(5$ = -inGSp . 



(10) 
(11) 



The Coriolis force is represented by the 2il x ^ term. 
SU is the hydrodynamic potential of the star, related to 
the Eulerian pressure perturbation 6p and gravitational 
potential perturbation by 



SU 



6p 
P 



(12) 



Assuming an e*"* time dependence, (IIOII and (|ll|l can 
be rewritten using the two-potential formalism of Ipser 
and Lindblom (this frequency convention differs from 
that used by Schenk et al jlQ] by a minus sign): 



1 



-4:TtGp^6{p) [SU + 5$] 



(13) 
(14) 



Note that for the rest of this paper the rotational 
frequency lo will be rescaled in terms of the parameter 
w = uj/2D, which takes on values between -1 and 1. (w is 
related to the k frequency parameter of Ipser and Lind- 
blom [III by u; = k/2.) 

Equation H14|) for the gravitational potential is 
Laplace's equation with a discontinuity on the boundary. 
The delta function describes the surface of the Maclaurin 
spheroid. Equation H13|l for SU is valid only within the 
spheroid. 

With two elegant coordinate transformations, Bryan 
transformed equations H13|) and H14|) so that they become 
separable and allow the boundary conditions to be im- 
posed. On the surface of the star these eigenfunctions 
correspond to waves that preserve the stellar volume. Im- 
posing the boundary conditions leads to dispersion rela- 
tions for the eigenfrequencies. 

The eigenfunctions found by Bryan ^3iE3 '^^^ 
pressed as Legendre functions in spheroidal coordinates 
(C, p) for the gravitational potential 



(5<I>i;|o 



a- 



a 



PuMe 



imcj) 



(15) 



(16) 



where the constant a is related to the eccentricity by 



t3 _ 



The hydrodynamic potential equation (|13(l is solved in 
terms of bispheroidal coordinates (^,/i): 



tn = 6v/(l-C2)(l-M2) 
(1 — ■w^)e'^ 



and 



P™(^o) 



(17) 



(18) 



where /3 is a constant related to a by the dispersion re- 
lations due to the boundary conditions. The constant 
Co = '^^Co /b'^d'^ is the coordinate value of the stellar sur- 
face in bispheroidal coordinates. Note that the Legendre 
functions are labeled using the indices n and m. This is 
to avoid confusion with the index / = n — 1 that is used to 
denote the spherical harmonics in the radiation reaction 
force, equation (j^IJ. 

The Lagrangian displacement ^ can be constructed 
from the hydrodynamic potential 13] equation (3.2). 



e 



1 - 

4^2' 



SUb 



where 



— 1 



(19) 



-V"/) (20) 



Finally the Lagrangian perturbation in pressure, equa- 
tion (15) of Ap = Sp + ^ ■ Vp can be written as 



Ap = p{SU + (5$) + ^ • Vp 



(21) 



using equation H12|l. The boundary conditions require 
that Ap = on the surface of the star. For the special 
case of the n = m + 1 r-modes it can be shown (see 
equation (7.11) of ^3) that Ap is identically zero over 
the entire star for any rotation rate. 



B. Properties of Linear Modes 

The coordinate transformations (|17|l used for the hy- 
drodynamic potentials depend on the mode frequency, 
w. Later, when the coupling coefficients are computed 
in Section IIVI it will be necessary to perform integrals 
involving three different modes; this requires a common 
coordinate system. It turns out that when transformed 
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back to cylindrical coordinates, the modes are nothing 
more than polynomials in w and z with frequency de- 
pendent coefficients, which makes the problem tractable. 
The construction of these polynomials is achieved by suc- 
cessive applications of recursion relations of the Legendre 
fmictions; this is presented in Appendix IXI The fact that 
the potentials can be represented in this way allows us 
to compute the large number of coupling coefficients be- 
tween the modes in a very efficient and accurate manner. 

Each eigenfunction can be labeled by three numbers 
(n, m, k) varying over the ranges specified below: 



n : 2 ^ oo 

m : ^ n — 1 

k : 1— >n — mifmT^O 
1 ^ n — m — 1 if m = 



Principal Legendre index 
Azimuthal index 
Frequency index. 



In order to estimate the number of modes for a given 
n, and the number of couplings among these modes, it is 
useful to introduce the index j which arranges the inertial 
modes consecutively: 



J 



The 



{n — l)n{n +1) (n — m — l)(n — m) 



6 



-k-l . (22) 



1 r-modes are 



j indices for the n — rr 
JR = (n - l)n{n + l)/6. 

To give some idea of the size and connectivity of our 
coupling network, observe that for a given maximum 
value of n = rimax , the number of included modes is 



N„ 



l)(n„ 



2) 



6 



6 



(23) 



The selection rules discussed in Section IVTl restrict the 
total number of couplings (K) to 



C. Frequency Spectrum 



(24) 



The frequency spectrum for the Bryan modes is set by 
the boundary conditions and is shown in Figure ^ 
The physical values of these frequencies are bounded be- 
tween — 2il and 2il. (Numerical values for the first few 
frequencies can be found in j^^l Table I ) 

For comparison, the oscillation frequency (/-mode) of 
an incompressible fluid sphere ( 12] equation (86)) is 



nGp 



8n(n- 1) 
3(2n+ 1) 



(25) 



More generally, the / and p modes of a compressible star 
have frequencies scaling as y/irGp. These frequencies are 
much larger than the r-mode frequencies even for mod- 
erately fast rotation. The large separation between these 
two frequency regimes makes the evolution of the r-mode 
difficult to follow in a hydrodynamical calculation. Any 



J , , — ^ V...- V. v.... 

■■*.■' r""* 



...-v V 



3 7 10 13 15 17 18 19 20 21 22 23 24 25 26 27 28 29 30 

n 

FIG. 1: Frequency spectrum for the inertial modes of an 
incompressible star, n — m + 1 r-modes denoted by The 
a;-axis ticks denotes the mode number n. The data points 
for the various m values follow each tick starting with the 
mode m = n — 1 nearest the tick mark and continuing with 
n — 2 ... 0. The j/-axis plots the frequency value w = u)/2Q,. 



direct fluid simulation requires a time step small enough 
to resolve the oscillations due to rapid pulsation modes, 
while studying the r-mode problem requires the evolu- 
tion of modes with much longer timescales. The large 
discrepancy between the frequencies of inertial and pul- 
sation modes makes it unlikely that the generalized r- 
modes would excite the pulsation modes. In all of our 
subsequent computations the pulsation modes will be ig- 
nored, allowing much larger time steps to be taken in 
mode-mode simulations, and enabling us to study the 
long term effects of mode-mode interactions. 



D. Normalization 

To facilitate the comparison between modes, all the 
linear modes are normalized by a factor ip such that the 
mode rotating-frame-energy at unit amplitude, e, has a 



MR'^n'^ where M = ^npR^ is 



fixed value of 2Eunit 
the mass of the star. 

Schenk et al. fl^ equations (K22) and (2.36) provide a 
general formula for the mode energy, valid at any rotation 
rate, namely 



1 2 2 

■ip W 



(fx pi* 



(26) 



These integrals can be computed analytically for an in- 
compressible star of any ellipticity. The computation in- 
volves a number of nonstandard integrals of Legendre 
functions [T3 | . In the spherical limit in units where R = 1 
the mode energy greatly simplifies to 



{n + m)! 



4172 (1 



'■) 2n + \ in — m)\ 



n(n+l)] . (27) 
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and the dimensionless normalization factor tp is 



1 — 2n + 1 {n — m)l 



n(n + 1) y 47r (n + m)! 



(28) 



The last square root factor is the standard normalization 
factor of Legendre functions to obtain spherical harmon- 
ics. 



For the case of the m = n — 1 r-mode the damping 
reaches its minimum value of 



^ = (2to + 3)(to-1) 
V 



(30) 



for a given n. Figure El displays the results for the shear 
damping. 



III. DISSIPATIVE EFFECTS 



B. Gravitational driving 



If the dissipative effects are small they can be modeled 
by adding a damping term jaca to the oscillator equa- 
tion (Q. For a detailed discussion of the validity of this 
approach see [l^ . 

The main damping mechanism in the higher order 
modes is viscosity, while gravitational radiation both 
damps and drives the low-order modes, with very little 
effect on the higher order modes. 



A. Shear viscous damping 

The damping rate, 7,, of a mode due to shear viscosity 
can be computed from the perturbed shear tensor |14| 



6a 



ah 



to be 



E„ 



2rj 



(fx Sa'^Hal^ 



where 77 is the shear viscosity coefficient. We use the 
coefficient proposed by [1^ for hot neutron star mat- 
ter, and explore temperature ranges varying from lO^K 
to lO^K: 



77 = 2x10^ 



P 



10^^ g cm 



9/4 /^o9 



T 



g cm "'^ s ""^ 



The geometric contribution {jn/v) of the individual 
modes is rather complex, having different scalings for a 
fixed n depending on the m value as well as the frequency. 



It varies from roughly 2n~ for the m 



1 r-modes to 



^n"^ for the m — modes. This picture is more compli- 
cated than the scalings given by Arras et al. in |0| based 
on the WKB analysis 

An analytic fit valid for all the modes motivated in 
part by the partial analytic results for the shear can 
be found for our numerically computed shear damping 
rates: 



1] 6 



[n + 3){n - 2) - 



2mw 



(29) 



Gravitational radiation strongly drives certain low- 
order modes, resulting in large scale fluid circulation 
while damping others. An accurate estimate of the driv- 
ing rates of gravitational radiation can be found using the 
procedure outlined by Ipser et al [l5l | and is summarized 
in the following equations: 



2w{2w - m)n'^ 
2E 



Y,Nii2w-mf n^' (I^AmP + \SJi„ 



l>2 



where SDim is the mass multipole (which vanishes in the 
spherical limit) 



6D, 



Im 



Sp r' 



(32) 



and 5Jim is the current multipole (which makes the 
largest contribution to driving in our model). 



5Jim. — — 



j r\2wnpi + 5pv) ■ y;f„* d^x 



(33) 



Yim and are the scalar and magnetic type vector 
spherical harmonics respectively. The quantity 



N, = 



AttG (/ + !)(; + 2) 

- 1)[(2/ + 1)!!]2 



(34) 



determines the strength of the coupling. The 0^'+-'^ factor 
in the denominator is the reason only large scale, small 
n modes such as the (n — 3,m ~ 2) r-mode are driven. 

The sign of the damping/driving rate is set by the 
2w(2w — m) factor. The eigenmode will be driven if this 
product is less than zero for that particular mode. This 
is the mathematical statement leading to the CFS insta- 
bility. 



C. Other sources of dissipation 

Bulk viscosity may be significant in neutron stars at 
high temperatures, particularly as a consequence of the 
appearance of hyperons at high densities 20] . The rate 
of dissipation of energy due to bulk viscosi ty J l4l| can be 
written using equations (114) and (117) of [lfl| as follows 
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FIG. 2: Damping due to shear viscosity of the inertial modes. The a;-axis ticks denote the mode number n. The data points for 
the various m values follow each tick, starting with the mode m = n — 1 nearest the tick mark and continuing with n — 2 . . . 0. 
The y-axis plots the geometric contribution of the shear 777/77. The inset displays the frequency and m dependence for n — 20. 
Lines join points with a fixed m value and are the calculated using the fit of equation 129|l . 



J (fx Quj^ \V-i\ 



Ap 





Ap 




P 



(35) 



Mathematically, this is zero for an incompressible star 



because —Eb oc F, 



0. We can adopt the attitude 



thatEq. 1(33) gives the dissipation rate to leading order in 
r|j~^, which we then regard as a parameter (which can be 
absorbed into the still-uncertain bulk viscosity coefficient 
C) . In the limit of slow rotation this integral simplifies to 



-Eb = 



p 



(36) 



where use has been made of equation H21|l and the fact 
that ^ scales like {l/fl^)VSU in equation ifTOI) . We have 
evaluated the integral (|36|l up to modes with n — 12 
and observe a strong frequency dependence in the results, 
modes with frequencies near zero being strongly damped. 
However, a deficiency of the incompressible star model 
for inertial mode perturbations is that Ap — for all 
n — m — 1 r-modes. 



Because we are mainly interested in the dynamical as- 
pects of three mode coupling in saturating the CFS in- 
stability of the r-mode, and since we cannot model bulk 
viscosity realistically or uniquely using an incompressible 
star model, we have chosen only to include shear viscos- 
ity in our network calculations. Viscosity plays two roles 
in our problem: (i) it is important in determining the 
parametric instability thresholds for interacting triplets 
of modes; (ii) it limits the range of modes that partic- 
ipate vigorously in the nonlinearly interacting network. 
Our goal during this investigation is to explore the char- 
acteristic behavior of mode-mode coupling in the simplest 
self-consistent model we can construct. The importance 
of bulk viscosity, as well as any other form of damping, 
can be examined retrospectively after we have explored 
the influence of damping of different strengths on our il- 
lustrative calculations for networks with shear viscosity 
alone. Accordingly, we reserve further comment on the 
omission of bulk viscosity to our next paper, where we 
report on the results of detailed numerical investigations. 
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D. Parameters of the Model 

All physical quantities in the calculation are made 
nondimensional by expressing mass in units of Mq , time 
in units of l/2f2, and length in units of the mean radius 
of the star, R. 

This choice of units allows the coupling coefficients 
to remain unchanged for all stars to first order in fi^. 
(Changing the rotation rate of the star influences the 
eigenfunction and this effect enters the coupling co- 
effiecients as second order in VP' .) It does however, mean 
that if the rotation rate of the star is changed both the 
viscous damping and GR driving rates are changed. The 
viscous damping rate changes because time is rescaled. 
The GR driving rate changes because enters nonlin- 
early into the formula for the GR driving. 

The dimensionless parameters in our simulations are 
chosen to correspond to the physical values of a star 
with i? = 10 km, M = IAMq and temperature rang- 
ing from lO^K to lO^K. The eigenf unctions considered 
in our calculations are those obtained in the limit of no 
rotation. However the gravitational driving is computed 
as if the star were rotating with a fixed rotation rate 
of 51 = O-iT^TipG corresponding to an eccentricity of 
e — 0.5. This rotation rate is fast corresponding to a 
physical frequency oi v — 695Hz or a period of 1.4 ms, 
which is appropriate since we are attempting to set an 
upper bound on the possible dynamics due to the r-mode 

I 



where ■~ab = CaAi^ Xabc = (^A.j^B.k^j.i and Ti = 
dhip/ d\n is the generalized adiabatic index govern- 
ing the perturbations. The notation Vs^abc} used for ex- 
ample in the term Sj^^V • ^^jj indicates that all cyclic 
permutations of the indices labeling the modes must be 
added, in other words V{abc} = Vabc + Vbca + Vcab ■ 

Note that we have changed the sign of $ from that used 
by Schenk et al f^, to be consistent with convention of 
Ipser and Lindblom and Chandrasekhar (17.] . 

We have chosen to investigate a particularly simple 
version of the nonlinear coupling, namely incompressible 
perturbations of a uniformly dense star. In principle, we 
could allow finite Fi for the perturbations, but we have 
made a tacit assumption that the perturbations and the 
star obey the same zero entropy equation of state with 
a single, large Fi. This is a convenient assumption for 
keeping our investigations of the nonlinear network rela- 
tively simple. 



driving. 

Our computed values for the damping and dri ving 
rates agree with those computed by Lockitch et al. [lj| 
et al for the case of a neutron star with R = 12.57 km, 
T = 1 X 10^ K, M = IAMq and = nGp, which are 
the values to which they scaled their results. 

The often-quoted growth time due to gravitational ra- 
diation of the order of 10'^ rotation periods is the maxi- 
mum possible growth rate that can be achieved if the star 
were rotating at breakup. Actual growth times could be 
orders of magnitude larger and increase like for the 
(n — 3,™ — 2) r-mode . 



IV. NONLINEAR COUPLING 

The linear growth of the unstable modes can be halted 
by various damping mechanisms. We wish to explore the 
damping due to couplings to other modes which even- 
tually transfer energy to small enough scales that it is 
efficiently dissipated. The pathways available for the en- 
ergy to dissipate are determined by the couplings and 
frequency detunings between the modes. We consider 
only the lowest order nonlinear three-mode couplings as 
computed by Schenk et al |10| . Their explicit expression 
for three-mode coupling as given by equation (4.20) fTd( 
is 



(37) 

I 

A. Nonlinear coupling for the incompressible star 

As the adiabatic index Fi —^ oo, the divergence of the 
Lagrangian displacement approaches zero (V-^ ~* 0). In 
the coupling coefficients, equation H37I) . these two terms 
appear as a product which has to be handled with care. 

To find an expression for p(Fi — 1)V'^ w pfiV-^ note 
that to first order the Lagrangian perturbations of den- 
sity (A/o) and pressure (Ap) obey Ap/p = — V-^ and 
Ap/p = TiAp/p (equations (114) and (117) of [l^) yield- 
ing 

pFiV-^ = ~Ap (38) 

Note that the term proportional to V • ^^V • $,b^ ■ 
vanishes because it is multiplied by only two powers of 
Fi. Thus the incompressible coupling coefficient reduces 



i^abc = ^ j d^x p{Ti-l)'E.{AB^ ■ ic}+P 



(ri - 1)^ 



d In p 



^■iA^- iB^ ■ ic 



\jd\p {xabc + Xacb) + ljd^xp ($|^eij<^^'^$c};y + ^\(Uc^■,^Jk 



8 



to 



K2 and K3 



Kabc 



i J (fx ^[AB^PC] - P iXABC + XACb) 

I J p[eiAeBS^''>'^c}■,^,+^\eBec'^■,vk 



(39) 



B. Coupling of Bryan modes 

The Bryan solution of the hnear modes is presented in 
terms of two scalar potentials 5$ equation l|14|l and 5U, 
equation from which the Lagrangian displacement 
^ is constructed in equation p9l) . It is in terms of these 
three quantities that we now construct the incompressible 
coupling coefficients (|39|) . 

Substituting the expression for the Lagrangian pres- 
sure perturbation into equation H39|) , the coupling coeffi- 
cient becomes 



KABC 



(fx p^^AB^Uc} 



~\ j '^^^ (^{AS^C}) -^P-P {XABC + XACb) 



(k5) 

(40) 



The domain of integration for ki, K2, and K4 is the vol- 
ume of the star, the Maclaurin spheroid of eccentricity e. 
The integrals for K3 and K5 are taken over all space; how- 
ever the factor p contains a step function that limits the 
integrals to the volume of the star and possibly a surface 
term. As written kabc has units of energy; computed 
values are appropriately non-dimensionalized by dividing 
by the mode energy, e. 

These terms can be rewritten in various ways that fa- 
cilitate, and check, the calculation of the coupling co- 
efficients. The only assumptions that are made during 
these manipulations are that the pressure p vanishes on 
the stellar surface, equation 0, and that the Lagrangian 
displacement has no divergence, V • ^ = 0. 

Schenk et al ^3 also give an alternative expression 
for the coupling coefficients in equation (4.27) of their 
paper. For completeness the correspondence between the 
coupling terms ki — kq and their expression is summarized 
in Appendix IdI 



The term K3 contains the second derivative of 5$. 
From equation (|14|l . this contains a delta function and 
yields a surface integral term of the same form as K5. 

Integrating (|14|l once gives rise to the boundary con- 
dition (equation (3.11) in [Tsf') 



47rGp2 
|Vp| 



{SU + (5$) 



STO 



^-[inGpn-^]^^, (41) 

where use was made of the Ap = boundary condition 
to obtain the second equation, S is a function that has a 
gradient equal to the outward directed unit normal vector 
Ua to the star and vanishes on the stellar surface. In 
particular 



P 



\Vp\ 



(42) 



is an example of one such function. 

The surface term in K3 can be computed by letting 

S<P = i/(-S)<5$ST0 + if (S)<5$EiO a (43) 

where H{x) is the Heaviside unit step function. Differ- 
entiating equation H43|) and using equation (|41|l gives 



(5$;y =i/(-S),5$ST0 ■^J + -ff (S)<5$SiO ■,^J 

- 4:TTGpS{J^)njnink£,'' 



(44) 



Using equation (|44|l . and integrating by parts once K2 
and K3 can be rewritten as an integral over the stellar 
surface only: 



k2 + k3 



d^A, ^TrGp^n.rik^^eBtc (45) 



As we show in Appendix El , the last term in equation 
is equal to (—3) x K5 



Term K4 is difficult to compute numerically as writ- 
ten. The two terms cancel to leading order in fi^ in the 
spherical limit, (e — > 0). Reorganizing the differentiation 
and using the Ap = boundary condition allows us to 
rewrite K4 as a complete divergence. The resulting sur- 
face integral is the same order in Q as the other terms in 
the coupling coefficient and can be integrated to yield a 
finite result. The details of the calculation are given in 
Appendix IbI with the result 

K4 = i y d'Ak p [{SUb + S<fB)ec-. + {6Uc + <5$c)Cb;.] e 

(46) 



^ / d^Ak ^{5UA + 5'^A)p'''CBihP■,^■,J 
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The result is asymmetric in the indices ABC, with no 
cychc permutation of indices imphed. (The expression 
can be symmetrized over the indices if desired by adding 
the two terms where the indices are cychc permutations 
of ABC and dividing the complete expression by three.) 



3. K5 

The term K5 scales as 0{il^) compared to the other 
terms that scale as O(il^). Appendix ICl shows that K5 
can be rewritten as the surface term. 

1^5^11 d^'x 47rGp2^(E) (^^^5^0"^'^^'^^) (47) 

In the spherical limit equation 147|) is identically zero. 

In the nonspherical case, this is the only term in the 
coupling coefhcients that cannot be calculated by means 
of our polynomial representation and term-by-term in- 
tegration (Section O. It does, however, yield an inte- 
gral that can be solved analytically as is shown in Ap- 
pendix El 



V. COMPUTATION OF COUPLING 
COEFFICIENTS AND OTHER INTEGRALS 

The computation of the coupling coefficients, shear 
damping, gravitational radiation driving or damping and 
normalization require the integration of quantities re- 
lated to 6U, (5$, ^ and their covariant derivatives to be 
computed over the volume of the background spheroid 
and its surface. The highly oscillatory nature of the func- 
tions for large n values precludes a direct numerical in- 
tegration. 

Our approach is to represent the functions analytically 
as far as possible. We store the linear potentials, SU, 
as polynomials in r and ro. These potentials can be con- 
structed directly from the recursion relations in z and vj 
as derived in Appendix ^ if the frequency of the eigen- 
mode is known. Integration of the coupling coefficient 
over the volume of the spheroid can then be achieved 
analytically. The coupling coefficient is computed as a 
polynomial and each term integrated individually. 

The analytic integral of a typical term in the coupling 
coefficient is obtained by working in the coordinates s 
and 9 which map the interior of the star onto the unit 
sphere and are defined by the coordinate transformation. 



zu = ReSsm9 z = RpS cos 9 (48) 

In these coordinates the volume element is 
di^x = R^RpS^ snv9 ds d9 d(j) and the integral over 



the star can be written as 

=Rl+''Rl+'' / e*(™-*+''"^+"^''^d</> 
Jo 

nTT p 1 

X / sin'^+i 9 cos'' 9d9 s^+^'+^'ds 
Jq Jo 

= 3_[^^_^_''^ '^(™^ + mB + mc) 

xB(i±l,i±i) + (49) 

where B{x, y) is the Beta function. 

All other integrals over the volume of the star, for ex- 
ample the normalization, damping etc. can be done in a 
similar way. The analytic result for normalization is used 
to check the accuracy of the method. 

In spite of using the analytic value for the individual 
terms in the integral, round-off errors become a serious 
problem. For the case where the coefficient is stored as 
a double precision floating point number, all significant 
digits in the normalization calculation have been lost by 
the time n > 16. Successive addition and subtraction of 
nearly equal terms cause the loss of precision. 

To overcome the loss in precision we made use of the 
GNU Multiple Precision Arithmetic Library GMP 4.1.2 
and worked to 32 digits accuracy. This resulted in the 
normalization integral differing from the analytic value 
by < 10^^'^ and coupling coefficients obeying their selec- 
tion rules listed in Table to within 10~^^ . The cost 
of this increased accuracy is speed: lO** couplings take 
approximately 4 days to compute on a Pentium 4, 2GIIz 
processor. 

VI. PROPERTIES OF COUPLING 
COEFFICIENTS 

The importance of a coupling in an oscillator network 
is determined by two factors: the magnitude of the cou- 
pling iKQ/j'yl and the detuning 

5Wal3-y = Wa — Wf} — (50) 

or frequency mismatch between the three coupled modes 
with frequencies Wq,, wp and w^. Modes with small de- 
tuning will ultimately play a far more important role in 
the dynamics of the problem than those with a larger 
frequency mismatch. Basically, a triplet with large de- 
tuning adds a very rapidly oscillating component to the 
nonlinear terms which on average will be zero. A small 
detuning results in a slowly varing contribution which 
may have a large accumulative effect and substantial in- 
fluence. 

The computed values for all the coupling coefficients up 
to a maximum value of n = 12 as a function of detuning 
are shown in Figure |21 
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-6t I I i I I 4 

-2-1012 
Detuning, 6 w 

FIG. 3: (Color online) The magnitude of the coupling coeffi- 
cient of a triplet of modes vs its detuning. In the background 
(red/light) are all the 766884 non-zero couplings for modes 
with n < 12. The foreground (blue/dark) shows the 3748 
direct couplings to the {n — 3, m = 2) r-mode. 

As mentioned in Schenk et al 10], the coupling coeffi- 
cients obey a set of selection rules. Our direct calculation 
verffies those stated in equations (5.11), (5.12) and (5.13) 
of their paper. An additional selection rule, the trian- 
gular selection rule on n, was discovered empirically. A 
complete list of selection rules obeyed by the coupling co- 
efficients in the spherical limit is summarized in Tabled 
All other coupling coefficients are found to be nonzero. 
We find that these selection rules persist for modes with 
nonzero eccentricity, at least for the few coefficients we 
checked. 



Selection Rule 


Statement of rule 


parity 


UA + ub + nc ~ odd number kasc ~ 


azimuthal, m 


niA + ruB + mc / ^ habc = 


triangular 


nc < |"-A — Jts ^ KAsc = 

UA + rlB < nc K.ABC = 



TABLE I: Selection rules obeyed by coupling coefficients 



Schenk et al [lOj also list a number of selection rules 
to 0(ri). If we consider ^ to be 0(1) then equation l|19|l 
implies that 6U and (5$ are 0(^1-^). Since all the coupling 
terms except k5 are composed of terms such as S,S,5U 
or ^^(5$, the lowest order coupling coefficient is 0{n^). 
Thus the additional selection rules in are trivially 
satisfied since they are only valid to 0(17). 

The empirically discovered triangular selection rule on 



n has also been seen in other calculations dealing with 2D 
planetary Rossby waves The coupling coefficients 

used in those calculations are less complicated than those 
considered in our case. 

In their paper on the saturation of the r-mode insta- 
bility p^ . Arras et al. make a momentum conservation 
assumption in their equation (A2) to obtain the scalings 
for a cascade solution. We find that all the couplings 
on which they based their calculation are zero since the 
momentum conservation assumption corresponds to the 
equality case of the triangular selection rule. However 
immediately adjacent to those modes are modes which 
give non-zero couplings and for which their WKB solu- 
tion may still be valid. 

We find that the contribution to the coupling coeffi- 
cients of terms containing perturbations to the gravita- 
tional potential are of the same order of magnitude as 
the hydrodynamic terms. This calls into question the 
applicability of the Cowling approximation, which is of- 
ten assumed to apply, for modes where the time scales 
are long and gravity has plenty of time to act. 



VII. DISCUSSION 

This paper provides the necessary analytic background 
and the computational approach required for calculating 
the couplings between the inertial modes of an incom- 
pressible perfect fluid star. These couplings determine 
the oscillator network of nonlinearly coupled modes that 
we shall use to explore the r-mode instability in neutron 
stars in Part II [Tl|. 

Although we have focused on the zero eccentricity limit 
in this paper, the techniques described are not limited to 
the spherical case, and remain valid for arbitrary rota- 
tion. An investigation into the effect of rotation on the 
eigenfunctions and frequencies, and resulting couplings 
and damping of modes will be an interesting study in its 
own right. 

A complete list of selection rules obeyed by the cou- 
plings is given in Tabled These rules along with the fre- 
quency spectrum shown in Figure ^ determine the basic 
nature and connectivity of the oscillator network. This 
structure will influence the nature of the dynamics of en- 
ergy transfer from the unstable large scale r-mode to the 
smaller damping scales and will set the maximum attain- 
able amplitude of the r-mode. The numerical integration 
of this network, and a detailed look at the non-linear dy- 
namics that result, will be discussed in detail in a Paper 
II . 11,] . The simplest dynamical system resulting from 
this network, the three mode Hamiltonian system, has 
already been used to analytically explain the dynamics 
of the catastropic decay of large amplitude r-modes ob- 
served in the hydrodynamical simulations It is ex- 
pected that investigation of the larger network will fur- 
ther increase our understanding of the evolution of the 
unstable r-mode. 

Thanks to Harald Pfeiffer for suggestions on the coding 
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of the classes and search algorithms used to manipulate 
and compute the coupling coefficients. 

This research is supported in part by NSF grants 
AST-0307273, PHY-9900672 and PHY-0312072 at Cor- 
nell University. 



APPENDIX A: GENERATING 
EIGENFUNCTIONS IN CYLINDRICAL 
COORDINATES 



Using equations lAip we derive the recurrence relations 
for the products of Legendre functions needed in equa- 
tions fn?)) and ifT^ . Define 



rm 

n 

rm 

ni 

jm 



[(2m- 1)!!] 



(A2) 



Our method of generating the potentials uses the re- 
currence relations of Legendre functions, equations (|A1|I 
to construct the potentials as polynomials in r and z. 



We shall also need the quantities 

G™ = xiP™(xi)p™i(x2)-fx2Pr(^2)pr-i(2^i) 



in-m)P:rix) = a;(2n-l)P™i(a;)-(n + m-l)P,7L2(a:) 
P;^'(x) = (-ir(2m-l)!!(l-.T2)f 



g: 



m+1 



p;„\i(x) = x(2™ + i)p™(a:) 



(Al) Using equations ljAl|l we can write 
_l 



(A3) 



{n - mfl-^ = XiX2[2n - l)2/™_i + {n + m~ 1)'C_2 " (2f^ " l)(^^ + m - l)G;r-i 



[n - m)Gl^ - {xl + xl){2n - l)C-i + 



[n + m — 1) 
(ri — 771 — 1) 



^ {{n + m - 2)G™ 2 - 2x^X2{2n ~ 2,)r^_^) 



(A4) 



Notice that the only combinations of coordinates 
that enter these relations for 5U , equation (|18(l with 
{xi = S,, X2 = fi) are 



^{l-xl){l-xl) 

{xl+xl) = l-^-^ 



J 



62 62^2 



XiX2 



bd 



(A5) 



The gravitational potential S^^io inside the star, equa- 
tion H15|l (xi = fi, X2 = iC) can be constructed using 



^{l~xl){l~xl) 



w 
a 



XiX2 



IZ 



(A6) 



APPENDIX B: TERM K4 AS A SURFACE TERM 

This appendix computes the surface integral version of 
K4 through successive integrations by parts and applica- 
tion of the boundary condition Ap = 0. The integrations 
by parts can be conducted by means of the manipulations 
listed in Table ITU Using Table HTl can be rewritten as 



follows. 

VXABC = £>! - £"2 + V5 -I- ^6 - ^3 

= Di - Da + P'4 + - Vii - Yxi - I^s + V^s + Vb 

(Bl) 

Let an over bar indicate that the B and C have been 
interchanged. Then 



VXACB = £>! - P'2 + V5 + H - 

= W^-D'2^Y^2+%~y'i . (B2) 

For the case of an incompressible fluid the pressure is a 
quadratic function in a;, y, and z and as a result Y\\ = 0. 
This allows K4 to be written as 



K4 



{Di -D2 + D^-D3 + Di- D2) (B3) 



which is a complete divergence and can be expressed as 
a surface integral over the boundary of the star. 

The pressure vanishes on the stellar surface, so, the 
terms Di, Di and can be ignored and 

K4 = i y d^x {-D2 + D4-D^) . (B4) 

Finally the Ap = boundary condition can be used to 
rewrite the nonzero terms: 
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Vi 
V2 

V3 
Vio 



Di - Vi 
V2 + V3 

D2-V4 
V^ + Ve 

D3 - V7 

Vs + Vg 

Di ^ Vio 
Vii + V12 



Vo 
Di 

V2 
V3 
D2 

Vs 

Di 

Vs 
V9 
Di 
Vii 



— PCA;j^B;k^C;i 
= PCA^B;k^C:i:j 



PCA^B;k^C;j 



= {p.rM^^Uh),, 
— P;r,k:it,At,B'iC 
= P;j;kCA^B:i^C 



TABLE 11: Terms generated when integrating pxabc by 
parts. Terms labeled D yield surface integrals while those 
labeled V are volume integrals. The left hand column shows 
successive steps in the integration process while the right hand 
column summarizes the terms. 



(B5) 
(B6) 



APPENDIX C: ANALYTIC COMPUTATION OF 



The computation of K5 requires a more detailed look 
at the background gravitational potential of a Maclaurin 
spheroid. We provide this background information in 
subsection IC II The integral expression for K5 in terms of 
the Lagrangian displacement is derived in subsection IC 21 
while the analytic evaluation of this integral in terms of 
Legendre functions is detailed in subsection IC 31 



1. Background - Ellipsoidal Figures of Equilibrium 

The background potential of a uniformly rotating star 
has been extensively studied. The results presented here 
are a brief summary of the parts of Chandrasekhar's book 
|l7j relevant to the computation of K5 

The gravitational potential obeying the equation 
— —AiiGp at an internal point Xi (cartesian coor- 
dinates) of a homogeneous ellipsoid with semi-axis is 



given by (equations (18), (40) of W^) 



where 



du 



A(af + u) 



(CI) 

(C2) 
u) and 



with A(u)2 = {a\ + u){al + u)(a§ 
/ = 010203 du/A. 

The gravitational potential at an external point can be 
found using equations (70), (71), (72) of and is 



^ x^ \ 
"Psio = 7rGpai0203 / -r 1 - 2^ — 



(C3) 



with A the ellipsoidal coordinate that is the positive root 
of the equation 



3 2 
xf 



E 



ai + X 



= 1 



(C4) 
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The gradients of the potentials are (no summation over 
repeated indices) 



'27rGpxiAi 
-2itG pxiaia2ci3 



du 



(C5) 
(C6) 



Note that on the surface of the star corresponding to 
A = the inner and outer potentials and their gradients 
are continuous. 



2. Integral Expression for Term K5 

The coupling coefficients require the term $;y7c multi- 
plied by the density of the star to be known. The poten- 
tial in the interior of the star is a quadratic function of 
x, y and z so it will vanish after taking three derivatives 
or <I>s|o -ijk — . The density outside the star also van- 
ishes. The only remaining contribution is due to a surface 
integral originating from the discontinuity in ^''-^ on the 
boundary which gives rise to a (5 function when taking 
the third derivative. To compute this term express the 
second derivative of the potential as 



Then 



(C7) 



(C8) 



where the first term is zero and the second does not con- 
tribute anything to the coupling coefficient. What re- 
mains is to calculate $Eio ^ *&sto -Aj on the boundary 
of the star: 



a2A(0) dxi 



A=0 

(C9) 



A=0 



Differentiating equation ljC4p with respect to and sub- 
sequently specializing to case of the stellar boundary 
gives 

3 9 



+ A dx^ ^ {a\ ^ \f 



dX 



A=0 



(CIO) 



Substituting ifUTUI) into (El gives 



= ATiGpn^n' 

Thus K5 due to the surface term is 
1 



(Cll) 



47rGp25(S) (ae'sCc^fe^.^ij) (C12) 



3. Analytic intergration of K5 

The integral for K5 (|C12|I can be computed analyti- 
cally by substituting the explicit expressions for ^ and 
5U and working in spheroidal coordinates on the surface 
of the star. After some manipulation equation ljC12|l 
reduces to: 



K5 



a2(4172)3 



/: 



V(l-'^') VCHCo) d^o 
• • - Jb [■■■]c 5{mA ■ 
1 



Co™ 
w 



TUB 



A 

mc) 



dp. 



(C13) 



This expression is symmetrical with respect to the mode 
indices, A, B and C and the term [• • -Js indicates that a 
similar expression as that contained in the square braces 
[• • -Ja for mode A should be included. Note that it is 
these three factors that determine the O(ri^) scaling of 
K5, while multiplication by the normalization factors, ip 
of each mode, cancels the il^ in the denominator of the 
constant factor. 

To compute the contribution of the integral, suppose 
u = V + w, < Z,n, m, and \u\ < /, \v\ < m, \w\ < 
n. Then the integral contribution rewritten in terms of 
partial fractions is 



dp ' ™ ' 



pu pv pw 

I 7n n 

2iCo 



1 



1 



p -iCo p + iCo 



(CM) 

The solution to the general case can be built up using 
induction. This method for calculatin g th e integral is 
based on the approach used by Gaunt |22| to calculate 
the integral of the product three Legendre functions. We 
first consider the two base cases: n = w and n = w + 1. 
Case n — w: Recall that 



p-?Ti _ / -1 \m ^ ™)- pm 

" " ^ ' {n + m)\ " 
(also holds for Q™) and use the integral [23l|p4 



(C15) 



dx 



(1 _ a;2yn/2 



P^" = (-2i)"(l - z2)™/2g™(z)z'' 



[z — a;) 

(C16) 

provided m < n; fc = 0, 1, . . . , n — m; z is in the complex 
plane with a cut along the interval (—1, 1) on the real 
axis. So 



pu pv pn 



{2ny. 



1 



1 



dn^ " 

1 (a^' + Co) 2»n!(2<o) J-i Vm - Ko M + Ko 

X [{l-p^)^Pl^] [(l-M^^-^p.] 

(C17) 
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Using Rodriques's Formula it can be seen that the second term is a polynomial of at 

most degree m — v. Soifm — v<l — u the integral can 
p / \ _ 1 d^jx ~ 1)" bg evaluated using equation (|C16|I to give 

~ 2"n! dx" 

P™(a;) = (-l)™(l-:.2)^^ (C18) 



1 pu pv pn 

df^ /2:rA = -(-2*)'gr(»Co)/'™(*Co)p„"(*Co)(i + (-i)'+™+") (ci9) 

1 [P- +^oJ 

Here use was made of = (-l)"P^'(z) and (3™(-2) = -(-l)"g;^'(z). 

Case w = n-l, M = 'y + n-l: Use P^^^ = (2n - l)xP^ri , so that 

pu pv pn—1 1 / 1 1 \ 



(2n- l)<o r / 1 1 



pM pn — 1 



(2zCo) J-iVM-^Co M + <0, 
= -(-2z)'Qr(^Co)P^(^Co)(2n - l)(^Co)P:^l\^Co) (l + (-1)'+'"+") 

= -(-2*)'Qr«o)p;„(*Co)pr'«o)(i + (-!)'+"+") (C2o) 

General case m — v < I ~ by induction: Suppose for all w < k < n 

"1 pu pv pw 

'^^'TY^ = -(-2z)™Qr«o)P;„Wo)/'r«o) (1 + (-1)'+"+'=) (C21) 
Then using the recursion relation, equation (|A1I) (a similar expression holds for Q™) it can be shown that 

' , PrPn.Pn ^ _J_ r l_\ pup. C^.H^.pu. _ n + W-1 

-lit^' + C^) (2zCo) li U - <o ^^ + ^cJ ' "U-^«^""' n-w ^ 



pu pv pw 



dll ' " = — — / — P^Pl ^lP^_^ ■ P, 

(2n - l)^Co / 1 , 1 

(n - w)(2<o) 7-1 \Ai - «Co M + <o 

+ f (-2^)'Qr«o)P^(^Co)P^-2«o) (l + (-l)'+™+«-2) " + ^ ~ V , 
\ n — w 

-(-2^)'Qr«o)P™(^Co) (l + (-!)'+'»+") f ^^""'^f P^-i - " + ^ ~ V , 

' \ (n — w) n — w 



n-2 



n-2 



-(-2z)'Qr«o)p;;(*Co)p„"'«o) (1 + (-1)'+™+") 



completing the proof by induction since the base cases of n y w and n — w + 1 were proved first. 



(C22) 



For the other cases m-v > l-u etc, there is nothing in APPENDIX D: (4.27) COUPLING 

this derivation that requires w, u to be positive. Thus you COEFFICIENTS OF SCHENK ET AL. 

could just relabel the Legendre functions appropriately 
and do the same calculation. The Q factor in the answer 
then always corresponds to the original P with largest 
difference between order and degree. 

Here we clarify the correspondence between the ki — K5 
of (|40() and those in equation (4.27) of which is re- 
produced in equation (|D1|) below. The terms containing 
V'^ have been omitted. The 5p term has also been re- 
placed using equation (|12|) . 
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KABC = - 2 y C{A^BKP^Uc})■,^J {W) 

-I J d'^X^{A^^SipS^C});rj («2) 

(Dl) 

In this case the domain of integration is all space and 
the pressure is taken to be p = S), that is, it van- 

ishes outside the star. All fields are taken to vanish at 
infinity, and complete divergence terms can be ignored. 
So doing integration by parts twice yields 7«r = ki and 
■^2 = K2, the step function in p restricting the integration 
to the volume of the star only. Note that 1^ = K3 and 
K5 = K5 trivially. 

The 7«4 term can be treated using the machinery de- 
veloped in Table ITU 

= - PXABC + Dl - D2 + Di - Dz 
+ V12 + 14 + ^9 (D2) 



and 

Vg-Vi2 = -pXACB+^l-^2+% (D3) 

so 

CA^BicP;i3k = ~p(XABC + Xacb) + Di - D2 + 

+01 -W2- D:i + V6 + VB+% 

(D4) 

All D terms can be ignored because the boundary is 
taken to be at infinity. Using the definitions of Table ^ 



14; + 14 + 1^6 = S^^sCc}?;- (D5) 

The only term that needs interpreting is p-^j 

p.,, = {H{~Y.)p),, = H{-Y.)p,, - 5{Y.)pn, (D6) 

but since p vanishes on the boundary, we are left with 
H{—Yj)p-j which restricts the integral to the volume of 
the star and yields an exact equivalence to the K4 term. 
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